!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!!
!! Simple Poisson problem for convergence tests.
!!
!! \n
!!
!! We consider the spatial domain \f$(-2\,,2)^{d_x}\f$ and the initial
!! condition
!! \f{displaymath}{
!!  f_0(\underline{x},\underline{v}) =
!!  \left(2\pi\right)^{-\frac{d_v}{2}}
!!  \exp\left(-\frac{1}{2}\sum_{\alpha = 1}^{d_v}
!!  \left(v^\alpha \right)^2 \right)
!!  \left( 1 + \frac{1}{4}d_x\pi^2\prod_{\alpha = 1}^{d_x}
!!   \cos\left( \frac{\pi}{2}x^\alpha \right) \right)
!! \f}
!! so that
!! \f{displaymath}{
!!  \rho_0(\underline{x}) = 1 + \frac{1}{4}d_x\pi^2\prod_{\alpha =
!!  1}^{d_x} \cos\left( \frac{\pi}{2}x^\alpha \right)
!! \f}
!! and
!! \f{displaymath}{
!!  \Phi_0(\underline{x}) = \prod_{\alpha = 1}^{d_x} \cos\left(
!!  \frac{\pi}{2}x^\alpha \right).
!! \f}
!! This can be used to test the Poisson solver at the initial time.
!<----------------------------------------------------------------------
module mod_poisson_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_poisson_test_constructor, &
   mod_poisson_test_destructor,  &
   mod_poisson_test_initialized, &
   test_name,        &
   test_description, &
   coeff_init,       &
   coeff_e,          &
   coeff_p

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter :: test_name = "poisson_test"
 character(len=*), parameter :: &
   test_description(1) = (/"Simple test for the Poisson solver"/)
 logical, protected :: mod_poisson_test_initialized = .false.

 ! private members
 real(wp), parameter :: pi = 3.1415926535897932384626433832795028_wp
 character(len=*), parameter :: &
   this_mod_name = 'mod_poisson_test'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_poisson_test_constructor(d)
  integer, intent(in) :: d(2)

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
          (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_poisson_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_poisson_test_initialized = .true.
 end subroutine mod_poisson_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_poisson_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_poisson_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_poisson_test_initialized = .false.
 end subroutine mod_poisson_test_destructor

!-----------------------------------------------------------------------

 pure function coeff_init(x,v) result(f0)
  real(wp), intent(in) :: x(:,:), v(:,:)
  real(wp) :: f0(size(x,2))

  integer :: l
  real(wp) :: ncoeff, tmp

  ncoeff = (2.0_wp*pi)**(-real(size(v,1),wp)/2.0_wp) ! normalization
  tmp = real(size(x,1),wp)*(pi**2)/4.0_wp
  do l=1,size(f0)
    f0(l) = exp(-0.5_wp*sum(v(:,l)**2)) * &
      ( 1.0_wp + tmp*product( cos(0.5_wp*pi*x(:,l)) ) )
  enddo
  f0 = ncoeff * f0 ! normalization

 end function coeff_init

!-----------------------------------------------------------------------

 !> Electric field, analytic solution
 !!
 !! The input argument \c x considers the fact that each component of
 !! the electric field could be defined on different nodes, so that
 !! the third index is used to differentiate among components. The
 !! result is thus
 !! \f{displaymath}{
 !!  E^\alpha(t,{\tt x(:,:,\alpha)}), \qquad {\tt x(:,l,\alpha)} =
 !!  \underline{x}_l^{E^\alpha}.
 !! \f}
 !!
 !! Notice that \f$t\f$ is ignored: this solution is exact only for
 !! the initial condition.
 pure function coeff_e(t,x) result(e)
  real(wp), intent(in) :: t, x(:,:,:)
  real(wp) :: e(size(x,1),size(x,2))

  integer :: id
  real(wp) :: tmp(size(x,1),size(x,2))
   
  do id=1,size(x,1) ! space dimension
    tmp       = cos(0.5_wp*pi*x(: ,:,id))
    tmp(id,:) = sin(0.5_wp*pi*x(id,:,id)) ! redefine the id-th comp.
    e(id,:) = -0.5_wp*pi*product(tmp,1)
  enddo

 end function coeff_e

!-----------------------------------------------------------------------

 !> Electric potential, analytic solution
 !!
 !! Notice that \f$t\f$ is ignored: this solution is exact only for
 !! the initial condition.
 pure function coeff_p(t,x) result(p)
  real(wp), intent(in) :: t, x(:,:)
  real(wp) :: p(size(x,2))

   p = product(cos(0.5_wp*pi*x),1)

 end function coeff_p

!-----------------------------------------------------------------------

end module mod_poisson_test

